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Abstract 

We consider the flow of a Newtonian fluid in a nano or microchannel with walls that have 
patterned variations in slip length. We formulate a set of equations to describe the effects on an 
incompressible Newtonian flow of small variations in slip, and solve these equations for slow flows. 
We test these equations using molecular dynamics simulations of flow between two walls which 
have patterned variations in wettability. Good qualitative agreement and a reasonable degree of 
quantitative agreement is found between the theory and the molecular dynamics simulations. The 
results of both analyses show that patterned wettability can be used to induce complex variations 
in flow. Finally we discuss the implications of our results for the design of microfluidic mixers using 
slip. 
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I. INTRODUCTION 



Several recent experiments 



□ □ □ □ 



report the measurement of large, shear-dependent 



liquid slip at partially wetting liquid-solid surfaces. While the origin of these dramatic 



violations of the no-slip boundary condition is still controversial 
to develop in how these effects may be exploited in microfludics 



interest is beginning 
. Microfludics is un- 



dergoing rapid growth with applications to chemical and biochemical synthesis |7[, and 
high-throughput synthesis and screening ^]. These applications require the manipulation of 
fluids in microchannels where flows are limited to very low Reynolds numbers. As a result, 
mixing in microffuidic devices tends to be diffusion dominated, requiring long channels and 
long retention times to achieve good mixing. As the scale of this technology continues to 
diminish the effects of low Reynolds numbers will become more significant. However, the 
effect of slip at channel walls also increases at small length scales so it is natural to ask 
whether the effects of slip can be used to overcome some of the disadvantages of laminar 
flow fj. 

To increase mixing rates it is necessary to induce transverse or circulating flows in a 
channel, increasing interfacial area between fluids or streamlines (for a recent review see 
This can be achieved by active mixers, which possess moving parts, but these can be difficult 
to fabricate and maintain. Passive mixers on the other hand achieve mixing by virtue of 
their topology alone and have no moving parts. Suggested designs for passive mixers include 



using channels with patterned topography \u}, 

in electro-osmotic flows jl^. Another possibility is to use chemically patterned channel 



111 ] , or channels with patterned surface charge 



151 ] have modeled the mixing of a 



surfaces. For example, Kuksenok and co-workers rJ M > 
binary AB fluid in channels patterned with A-like and B-like regions. 

Yet another approach might be to use patterned wettability to induce variations in slip. 
Slip is often characterised by a slip-length 5, which is the distance at which the fluid ve- 
locity at a surface (i.e. the slip velocity) vanishes if it is linearly extrapolated beyond the 
surface. Measurements of slip lengths do vary widely but some groups have reported slip 
lengths of several microns l|. It is common to invoke the formation of nanobubbles at 
the hydrophobic surface explain such large slip lengths. However there is still much 

disagreement about the magnitude of slip that can be induced although lengths of tens to 



hundreds of nanometers seem to be more typical 



Furthermore other factors such as 
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surface roughness and surface contamination do seem influence slip length measurements 
considerably jj, [3]. Lauga and Stone [itJ have recently considered the effect of patterned 
no-slip and no-shear stress regions in pressure-driven Stokes flow in a cylindrical geometry 
where the no-shear stress regions model the presence of micro or nanobubbles on channel 
walls. From this they derive an effective macroscopic slip-length which indeed is found to 
depend on shear-rate and geometry. 

Here we will assume that the variations in wettability can produce variations in the slip 
of the flow at the channel walls. Molecular dynamics simulations of flow past hydrophobic 
surfaces [isl, Q, 13 certainly demonstrate a strong relationship between wettability and slip, 
although the slip lengths found tend to be of the order of a few of tens molecular diameters 
at most. However, as noted above, the formation of nanobubbles at hydrophobic surfaces 
may well be able to induce very large slip lengths: thus large variations in wettability on a 
surface might be expected to produce large variations in slip length. 

We begin by studying a Newtonian flow in a simple channel with a slip boundary condition 
characterised by a slip-length 5 that varies in space i.e. 5 = 5(x). In the first instance we are 
interested in seeing the effect of a variable slip length on the flow, and secondly, in evaluating 
the possible exploitation of such effects in designing a microfluidic mixing device. As this is 
a preliminary study, we will approach the problem analytically, using a pertubative scheme 
to satisfy the slip boundary condition. 

In section 4, we will use molecular dynamics simulations of the flow of a Lennard- Jones 
fluid between two plates. The interaction between the plates and the fluid will be allowed to 
vary in space in order to test the predictions of the analysis in sections 2-3. We conclude with 
a discussion of the implications of our findings here for the design of mixers in microchannels 
using chemical patterning. 

II. EQUATIONS FOR FLOWS WITH SPATIALLY VARYING SLIP 

We start with the Navier-Stokes equations for a viscous incompressible fluid: 




(1) 



V-u = 0, 



(2) 



3 



where u is the velocity field, p the fluid pressure, p is the fluid density and /x is the fluid 
viscosity. 

We consider a pressure driven flow in a two-dimensional channel geometry corresponding 
to flow between two plates as illustrated in figure The channel has length L and width 



2w. At the channel walls we have Navier slip boundary condition 



2l|: 



r)ii 

u(±w) = T$-f(±w) (3) 

where u is the longtunidal velocity component (x-direction). The transverse velocity com- 
ponent v (y-direction) satisifies v(±w) = at the walls. At the channel exit and entry we 
prescribe the pressure to be p and Pl respectively giving a pressure head across the channel 
of Ap = po -p L . 

The solution to (HHHJ) is 

u = ^U + 2w5-y 2 )= u(l + 2--£) (4) 
v = (5) 

P = p in - Ap (^j (6) 

where U = ^ is the maximum fluid velocity in the absence of slip. 
We will now allow the slip length to vary in the x direction i.e. 

u(±w) = tS(x) — (±w). (7) 

oy 

Specifically, we will consider the following slip boundary condition: 

t (±w) = T 5 (l + ae ikx ) ^{±w). (8) 



u 



If a <C 1 then we can apply a perturbative approach 

u = uo + aui + O (a 2 ^ (9) 
v — Vq + avi + ... (10) 
p = po + ap 1 + ... (11) 

where (wo,fo,Po) solve the constant slip-length boundary condition problem (JT]|HJ). The 

equations for the first-order corrections in a are then given by 

/ dui du \ dpi 2 
P \u ^- + vi— )= -— +yUV Ui (12) 
\ ox oy J ox 

P^ = ~^+^ Vl (13) 



and 



dUl + ^ = (14) 



dx dy 
with boundary condition 

Ul (± W ) = 2 -^e^TS d -^(±w) + 0(a). (15) 
fxL dy 

The boundary condition immediately suggests the solution ansatz U\ = e lkx f{y). Inserting 
this into equation (j!4jl . we find that 

Vl = -ik e ikx h{y) (16) 

where h'(y) = f(y) and h(0) = since fi(0) = by symmetry. 

We can now eliminate pi from (J12|) and (|13|) to obtain the following ordinary differential 
equation for h(y): 

-0+(^+^0-^+H^ + f))^° (17) 

where v = fi/p is the specific viscosity. In terms of h the boundary condition (fTHjl becomes 

^(±w) = 2U- =F 5p,{±w) + 0(a). (18) 
ay w dy z 

We note that differential equation (fTTj) is homogeneous, so the magnitude of h will be set 
by the boundary condition (fTH|) . Further when a C 1, (fTTjl and (fTK|) form a quasilinear 
boundary value problem. In the following section we will examine the solution to this 
problem in a number of limiting cases. 



III. RELEVANCE TO MICROFLUIDIC DEVICES 

At this stage we will introduce some scales into the problem. As the effects of boundary 
slip on the flow scale as 5/w (see equation HJ), at widths substantially greater than the slip 
length, a surface with patterned wettability will have an insignificant effect on the flow. 
With values of the slip-length of up to several microns reported in the literature |jj , we will 
confine ourselves to discussion of channels with widths w < 10 pm. 

Obviously the wavelength of the slip variations 2ir/k is bounded above by the length of 
the channel L. The wavelength is bounded below by the minimum size on which the slip 
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can be patterned. While microcontact printing |22j or photolithography of hydrophilic or 
hydrophobic self-assembled monolayers might be limited to wavelengths greater than several 
hundred jum, in principle, it is still of interest to consider the limit as 2ir/k ~ 10 nm. This 
might achieved using a combination of self-assembly by block copolymers and lithography for 
example |23|. Hence it is reasonable to consider patternings that satisfy 10 7 nm > 2n/k > 10 
nm. 



A. Slow flows with fine patterning 

In the limit where -C 1 and <C k 2 w 2 then equation (fTTj) reduces to 

= (i9) 

For instance, in the case of water which is flowing at 0.01 ms" 1 in a 10 /im-width channel (i.e. 
1 nanoliter per second) then U/kv -C w 2 k 2 <C 1 for 1/k <C 20/im. Note that equation ()19j) 
is real (whereas the equation (fTTj) is complex) so the variation in longitudinal flow velocity 
is in phase with the variations in slip length while the variations in transverse flow velocity 
are 90° out of phase with the variations in slip length (recall equation (fTo^ ). 
The solution to ()19|) with boundary condition ()18j) to order a is given by 

h( ) = jj ( -\ ( w ~ y) sinh k(w + y)- (w + y) sinh k{w -y) 
[V) ~ \w) smh 2kw + 2k5 cosh 2kw - Ak(w + 5) [ ' 

where we recall that v\ = —ike lhx h(y) and u\ = e lkx h'(y). Figure El shows h(y) for kw = 10, 

1 and 0.1. It is clear from ()19)) and figure El that the magnitude of kw controls the variation 

away from the walls in h(y) and hence in v\ and %. With kw <C 1, then the transverse 

velocity induced V\ is confined to very near the walls. Indeed, from figure El we see that to 

maximize both the magnitude of h(y), and its penetration towards the centre of the channel 

we should choose kw ~ 1. Similarly, figure 01 shows h'(y) for kw = 10, 1 and 0.1. 

Figure |U shows a flooded contour plot of the variations in both velocity components (u\ 

and v±) along a long channel (L = 20w) with kw = 1, and figure El shows a vector plot of the 

velocity components in a shorter channel (L = nw) with kw = 1. Note that where the slip 

at the channel wall is high, the flow velocity increases at the channel walls, but decreases in 

the center of the channel. Likewise, where the slip is low, the flow velocity decreases at the 
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channel wall but increases in the center of the channel. Between the peaks and troughs in 
slip, transverse flow is generated away from or towards the channel walls. 

We can also look at square-wave variations in slip length, utilizing the Fourier series for 
a square wave of wavelength A = 2ix/k: 

/(x)= 4g sin((2„ pi) 

71 n=l ZTl 1 

Since equation is linear we can solve for each Fourier mode and resum to obtain the 
solution for a square wave variation in 5(x). Doing so gives 

4 00 

vi = -Y / kh(k,y)cos((2n-l)kx). (23) 

n n=l 

Figure |H1 shows a flooded contour plot of the velocity components v\ and u\ for a square 
wave variation in 5 with kw = 1. 



B. Slow flows with larger scale patterning 

Now we consider the situation where 1 ^> ~ k 2 w 2 . For instance, in the case of water 
which is flowing at 0.01 ms" 1 in a 10 /zm-width channel (i.e. 1 nanoliter per second) then 
U/ku ~ w 2 k 2 <C 1 for 1/k ~ 20/im. This corresponds to a slow flow with spatial variations 
in slip length occuring on scales greater than the channel width. Now equation (|T7jl reduces 
to: 

d 4 h „,od 2 h / „ . U \ , 
The solution to this equation with boundary conditions ()18j) is 



h(y) 



U (^^j (sinh X + y sinh A_w — sinh X^y sinh A^ 



w) 



A + cosh A+u; sinh \-W — A_ cosh A_w sinh \_w + 5(A+ — A?_) sinh A+w sinh A_w 

(25) 



where 



^t-'Jk- < 26) 

Note that when w F k -s u —> 0, the expression (|23j) for /i(y) reduces to the expression (j2~Uj) from 
the previous section. In fact it is instructive (although tedious) to write (J23|) in the form of 



([20)1 plus corrections in w F k 3 u - Doing so we can write h (y) as: 

ft fr>=* w fr> + ^ kM fr> +0 ((dk)l (27) 

where h^ 1 ' (y) is given by equation EDI 

(y) = -—- (Qk(w 2 — y 2 ) sinh ky sinh kw + (w + y)(k 2 (w + y) 2 — 3) sinh k(w — y) 
12k v 

— {w — y){k 2 {w — y) 2 + 3) sinh k{w + y) 

(y) (8k 3 w 2 (35 + w) + 6k 2 cosh 2kw + 3(A5k 2 w - 1) sinh 2kw)) (28) 

and 

k = sinh 2kw + 2k5 cosh 2kw - Ak(w + 5). (29) 

Note that the first order correction in aj^ is purely imaginary. Thus it introduces a phase 
lag the response of the fluid to the slip at the walls (moving it downstream) and increases 
the magnitude of h(y). This is shown in figure [7| which compares the velocity u\ in the 
centre of the channel for A5 = 5sinkx (kw = 0.1) for w F k a u — and w H k z v = 15 (a large 
value of R 3 is chosen here so that the effect of this term is easily visible). 



IV. MOLECULAR DYNAMICS SIMULATIONS 



To study the effect of spatially varying wettability on flow in a channel at a molecular level, 
we have followed the approach of Barrat and Bocquet We consider a Lennard- Jones 

fluid with atomic mass m confined between two walls with periodic boundary conditions 
applied in the plane of the walls. The walls consist of fixed Lennard- Jones atoms and 
interact with the fluid via a modified Lennard- Jones potential of the form: 



4e 



12 



a 

r. 



a 



C fs 



(30) 



where < Cf s < 1 controls the degree of wettability of the walls [18[. Note that the fluid 
atoms also interact according to potential (J30)) with Cff = 1.2. Here we will consider flows 
where Cf s = Cf s (x) to model the effect of chemical patterning of the channel walls. 

We used a simulation cell containing 6750 fluid atoms within a volume of approximately 
(20a) 3 . The temperature was controlled using Nose-Hoover thermostat 24\ on the velocity 
component of the fluid atoms parallel to the channel walls but perpendicular to the imposed 
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flow direction (in figure ^ this is the direction into the page). Flow can be induced by 
applying a body force to the fluid atoms in a direction parallel to walls giving a Poiseuille- 
type flow, or by dragging one of the walls past the fluid which induces a Couette flow. Here 
we only consider the former as our intention is to make a comparison with the pressure-driven 
flows of the previous section. 

When Cf s = 1.0 everywhere the flows are well- approximated by solutions to the incom- 
pressible Navier-Stokes equations (012)) with a no-slip boundary condition, although density 
variations occur near the walls due to the well-known tendancy for fluid atoms to layer at 
a solid interface. Furthermore, when 0.5 < c/ s < 1.0 but is constant everywhere, we find 
that the flow is reasonably well-approximated by solutions to the Navier-Stokes equations 
with a slipboundary condition (j3J). Our simulations are in good agreement with Barrat and 
Bocquet 

We now consider the simulation of flows in a channel with patterned slip length. The 
equation for Cf s on the channel walls is: 



c fs 



0.9 ifsinfcx>0 

(31) 

0.5 if sin kx < 



where kw = 7r (i.e. the wavelength of the pattern is 20 a, which is the width of channel). 
Note that the mean value of Cf s is 0.7. Our simulations show that such a patterning does 
indeed induce a variation in slip length along the channel walls. For instance, as illustrated in 
figure IHlfor a simulated flow with peak flow longitudinal flow velocity U = 1.30 (e/m) 1 / 2 , we 
calculated an effective slip length of 5 = 9.1 a by fitting a parabolic profile U(l+S/w— y 2 /w 2 ) 
to the longitudinal velocity profile in the solvophilic region (x > i.e. where Cf s = 0.9). In 
the the solvophobic region (x < Cf s = 0.5) we calculated an effective slip length 5 = 13.0 a. 
Similarly, for a simulated flow with peak flow longitudinal flow velocity U = 0.4 (e/m) 1 / 2 , 
we calculated an effective slip length of 5 = 3.6 o in the solvophilic region [x > i.e. where 
Cf s = 0.9). Likewise in the the solvophobic region (x < c/ s = 0.5) we calculated an 
effective slip length 6 = 6.7a. 

Furthermore, these variations in effective slip length induce transverse flows as anticipated 
in the previous sections. Figure 01 shows the time-averaged transverse velocity of a flow, 
with peak longitudinal flow velocity U = 0.4 (e/m) 1 / 2 . The peak transverse velocity is 
V = 0.03 (e/m) 1 / 2 . Re gions with light shading indicate flow in the y-direction and regions 
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with dark shading indicate flow in the negative ^/-direction. Note that the variations in 
v(x, y) are 90° out of phase with the variations in Cf s as predicted by our analysis in section 
3. To compare the magnitude of the variations in v from the simulation to the theory of 
the previous sections, we use the effective slip lengths calculated above. Thus 5 = 5.2 a 
and a = 0.3 in equation (jSJ). Figure ITU1 compares the theoretically expected variation in v 
at x = across the channel for a square wave variation in slip length (see equation (|23p) 
to the time-averaged simulated variations. It is seen from the comparison that the theory 
underestimates the peak values of v by a factor of 2-3. 

Figure HP shows a faster flow with peak longitudinal flow velocity is U = 1.30 (e/m) 1 / 2 
and peak transverse velocity is V = 0.060 (e/m) 1 / 2 . Regions with light shading indicate flow 
in the ^/-direction and regions with dark shading indicate flow in the negative ^/-direction. 
Note the downstream phase lag in the variations in v(x,y) with respect to the variations in 
v(x, y) in the slower flow shown in figure H3 We have not made a direct comparison of this 
phase lag with the predicted phase lag in equation (|2*7jl as we were unable to solve the full 
equation for h (fTTj) for fast flows analytically. However, once again we find that the theory 
underestimates the peak values of v by a factor of 2. 



V. DISCUSSION 



The molecular dynamics simulations in section 4 demonstrate that patterned wettability 
will induce patterned variations in slip length. While this is not surprising given the demon- 
strated link between wettability and slip in other molecular dynamics simulations 
supports the use of the boundary condition (JJJ) in evaluating the effect of patterned wetta- 
bility on flow. Furthermore we found a strong qualitative agreement between the molecular 
dynamics simulations and the approximate analytic solutions developed in sections 2-3, al- 
though the theory tended to underestimate the magnitude of the variations in flow due the 
patterning by a factor of 2. This underestimation may in part be due to the way the theory 
was fitted to the simulations (i.e. by fitting effective slip lengths in the different channel 
regions). However, the theory also assumes the fluid is incompressible, whereas substantial 
variations in fluid density can occur at the walls. In particular, a reduction in the density 
of the fluid near the solvophobic region of the wall relative to the solvophilic region of the 
wall, as observed in the MD simulations, would tend to enhance the transverse variations 
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in flow. Finally, we note that in our analysis in sections 2-3 we only solved the slip bound- 
ary condition to order a (the relative variation in slip length). In the molecular dynamics 
simulations conducted here a was found to be 0.2-0.3. 

In any case our calculations here have demonstrated that in an incompressible Newtonian 
fluid, changes in slip length can induce flow transverse to the walls in a nano or microfluidic 
channel. Further our calculations suggest that these transverse flows are maximised if the 
patterning of slip takes place on a wavelength A ~ w. Thus it certainly appears that 
patterned slip could be used to induce mixing in the same way as patterned topography 
(e.g. such as the asymmetric herringbone pattern studied in Ref jlQj]). Figure IT2l suggests 
some possible patternings that could be used for mixing. However, we note that slip also 
changes the velocity profile in a channel (e.g. see figure |HI). These changes in profile will no 
doubt alter the effect of dispersion on mixing in a channel. While our perturbative approach 
is not valid in the limit of large variations in slip-length which would be most effective for 
mixing, we would expect the flows to be qualitatively similar. Cottin-Bizonne et al j^] have 
calculated effective slip-lengths in a half-plane geometry for flows over no-slip and partial 
or full slip patterned regions, although they have not examined how this alters the velocity 
profiles. 

We also note that surfaces with switchable wettability have recently been demonstrated 
2(| . This switchability suggests the interesting prospect of a slip length which is time and 
space dependent i.e. 6 = 5(x,t). The approach outlined in section 2 can easily be adapted 
to consider this situation. If we imagine a traveling wave variation in slip length 5(ut + kx), 
then in a frame comoving with this wave, the effects on the flow will appear similar to those 
of speeding up (or slowing down) the flow. Thus the response of the fluid to a rapidly 
changing time-dependent slip will lag these changes in slip (as the fluid response does for 
the fast moving flow in figure |7|). We will consider this problem in more detail in further 
work. 

VI. CONCLUSION 

We have considered the flow of a Newtonian fluid in a channel with spatially varying 
surface properties. Using a pertubative approach we derived equations that describe flow 
in a channel with patterned variations in slip length. We also examined flows in a channel 
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with varying wettability using molecular dynamics. The simulations demonstrated that the 
variations in wettability induce variations in slip. Good qualitative agreement was found 
between the molecular dynamics simulations and the approximate analysis of the Navier- 
Stokes equations. 
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FIG. 1: The two-dimensional channel geometry. 



14 




FIG. 2: The function kh(y) ~ v\ is shown in a channel for values of kw=0.1,l,10 respectively. We 
have taken 5/w = 1. 
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FIG. 3: The function h'(y) ~ u\ is shown in a channel for values of kw=0.1,l,10 respectively. We 
have taken 5/w = 1. 
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FIG. 4: Contour plot showing A5(x)/5 = cos(kx) and the corresponding variations in v\(x, y) and 
ui(x,y) in a channel for kw=l. Regions with dark shading indicate negative velocity and regions 
with light shading indicate positive velocities. We have taken 5/w = 1 and the channel length is 
L = 20w. 
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A5(x)/5 




x/w 



FIG. 5: Vector plot showing A5(x)/5 = cos(/cx) and the corresponding velocity vector 
(ui(x,y),vi(x,y)) in a channel for kw=l. We have taken 6/w = 1 and the channel length is 

L = TTW. 
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FIG. 6: Contour plot showing a square wave AS(x)/S and the corresponding variations in vi(x, y) 
and ui(x,y) in a channel for kw=l. Regions with dark shading indicate negative velocity and 
regions with light shading indicate positive velocities. We have taken S/w = 1 and the channel 
length is L = 20w. 
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FIG. 7: Plot of the change in velocity down the channel for w $ k3u = 15 (solid line) and for 
^ = (dashed line) for kw = 0.1. We have taken 5/w = 1. The effect of this term is to 
cause a phase lag in the velocity corrections away from the walls (effectively shifting these changes 
downstream) and to increase the magnitude of these corrections. 
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FIG. 9: A plot showing the square wave Cf s (x) boundary condition with kw = it imposed on 
the walls of the molecular dynamics simulation, the effective slip lengths 5 e Q induced by Cf s and 
a corresonding contour plot showing the variations in v(x,y) = v\(x,y). The channel width is 
2w = 20<7 with periodic boundary conditions applied at x = ±10<r. The peak longitudinal flow 
velocity is U = 0.4(e/m) 1 / 2 and the peak transverse velocity is V = 0.03(e/m) 1//2 . Regions with 
light shading indicate flow in the y-direction and regions with light shading indicate flow in the 
negative y-direction. Note that the variations in v(x,y) are 90° out of phase with the variations in 

Cfs- 
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FIG. 10: A comparison of the transverse velocity v at x = across the channel from the simulation 
in figure ED and from theory. The theory has been fitted to the simulation data by calculating an 
effective slip length 5 = 6.7 a across the solvophilic region (x > i.e. where c/ s = 0.9) and an 
effective slip length 5 = 3.6 a across the solvophobic region (x < cj s = 0.5). Thus 5 = 5.2a and 
a = 0.3. It is seen from the comparison that the theory underestimates the peak values of v by a 
factor of 2-3. 
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FIG. 11: A plot showing the square wave Cf s (x) boundary condition with kw = ir imposed on the 
walls of the molecular dynamics simulation, the effective slip lengths 5 e Q induced by Cf s and the 
corresponding variations in v(x,y) = vi(x,y) in a channel. The channel width is 2w = 20a with 
periodic boundary conditions applied at x = ±10c The peak flow velocity is U = 1.30(e/m) 1//2 
and the peak transverse velocity is V = 0.060(e/m) 1 / 2 . Regions with light shading indicate flow in 
the y-direction and regions with light shading indicate flow in the negative y-direction. Note the 
downstream phase shift in the variations in v(x,y) with respect to the variations in c/ s , especially 
in comparion with figure El 
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FIG. 12: Suggested designs for mixing devices. The light regions would be coated in such a way 
as to induce a large slip length (say with a super hydrophobic coating), while the dark regions 
would be coated to induce a small slip length or no slip (say with a superhydrophilic coating). 
More complicated patterns may enhance the mixing, provided the patterns are on a length scale 
comparable to the channel width. 
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